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ABSTRACT 

We develop a lattice Boltzmann equation method for simulating multi-phase im- 
miscible fluid flows with variation of density and viscosity, based on the model pro- 
posed by Gunstensen et a/|| for two-component immiscible fluids. The numerical 
measurements of surface tension and viscosity agree well with theoretical predictions. 
Several basic numerical tests, including spinodal decomposition, two-phase fluid flows 
in two-dimensional channels and two-phase viscous fingering, are shown in agreement 
of experiments and analytical solutions. 
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1 Introduction 



It has recently been shown that lattice gas computational fluid dynamics, including 
lattice gas automata (LGA) [Ij] and lattice Boltzmann (LB) equation methods, pro- 
vide alternative numerical techniques for solving the Navier-Stokes equations, multi- 
phase fluid flows |, |5j and a variety of other fluid systems || [7|. The parallel nature 
of these newly developed schemes, adapted from cellular automata, affords an easy 
implementation of fast, efficient and accurate simulations on parallel machines. 

To solve the incompressible fluid equations by traditional numerical methods such 
as finite differences or finite elements, one must deal with a Poisson equation for 
the pressure term that is induced by the continuum condition and the momentum 
equation. This fact has brought to light a crucial difference between lattice gas 
methods and traditional algorithms, since a lattice gas solver is only required to 
solve the kinetic equation of the particle distribution. All other quantities, including 
density, velocity and energy can be obtained by a macroscopic averaging process 
through the distribution function. Therefore, the pressure effects on the momentum 
equation are controlled by an equation of state. 

Strictly speaking, lattice gas automata should be restricted to solving compressible 
fluid flows. The incompressible limit theory of lattice gas methods || and their numer- 
ical simulations^, however, predict and demonstrate that such a lattice gas solution 
at low Mach number produces results comparable with those of the incompressible 
Navier-Stokes equations for a wide variety of problems and with very small numerical 
errors. Also, two-phase surface dynamic boundaries and wall boundaries are far more 
easily implemented with a lattice gas method compared to direct simulations of the 
incompressible Navier-Stokes equations |3|, (F|. 

Rothman and Keller || were the first to extend the single-phase lattice gas model 
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proposed by Frisch, Hasslacher and Pomeau|l[] to simulate multi-phase fluid prob- 
lems. Colored particles were introduced to distinguish between phases, and a nearest- 
neighbor particle interaction was used to facilitate interfacial dynamics, such as 
Laplace's formula for surface tension. Later, Somers and RemQ], and Chen et a/|| 
extended the original colored particle scheme by introducing colored holes. It has 
been shown that the colored- hole lattice gas method extends the original nearest 
neighbor particle interaction to several lattice lengths, leading to a Ukowa potential. 
Moreover, the colored-hole scheme carries purely local information in it's particle col- 
lision step, reducing the size of the look-up table in the algorithm and consequently 
speeding up the simulation. 

Although two-phase lattice gas algorithms are able to produce interesting surface 
phenomena, they are difficult to compare quantitatively with experiments and other 
numerical simulations due to their noisy characteristics induced by particle fluctua- 
tions. The lattice Boltzmann model proposed by McNamara and Zanetti0, however, 
solves the kinetic equation for the particle distribution instead of tracking each par- 
ticle's motion, as is done in lattice gas methods. Nonetheless, lattice gas algorithms 
will be superior to lattice Boltzmann models for simulating underlying physics and 
modeling microscopic dynamics such as correlation effects and phase transitions, since 
they contain more information about the microscopic behavior of particles. On the 
other hand, to provide a numerical method for solving partial differential equations 
governing macroscopic behavior, lattice Boltzmann schemes will be at least as good 
as lattice gas models. The finite-difference nature of lattice Boltzmann methods not 
only can simulate the macroscopic equations more efficiently and accurately, but also 
preserve some of the advantages of lattice gas models, such as their parallel computing 
nature and their ease of boundary implementation. 

Combining the lattice Boltzmann model of McNamara and ZanettiH and the 
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two-phase lattice gas model of Rothman et al, Gunstensen et al proposed a lat- 
tice Boltzmann method for solving two-phase fluid flows. An important contribution 
of this model is the introduction of a perturbation step (shown in detail below) so 
that Laplace's formula at an interface can be approximately recoved. This lattice 
Boltzmann method has been used in several applications fT0|| , however it has a few 
fundamental problems. First, the model does not solve the exact two-phase fluid 
equations - although Galilean invariance is recovered by the proper assignment of 
rest particles, the equation of state remains to be velocity dependent [0 . Secondly, 
the model uses a fully linearized collision operator, which involves a 24*24 matrix 
multiplication at each time step and position in three dimensional space, reducing 
computational efficiency Finally, the model is restricted to two-phase fluid flows 
having the same densities and viscosities. In this paper, we extend their model by 
using a single-time relaxation lattice Boltzmann model. With the proper choice of the 
particle equilibrium distribution function, we are able to recover the incompressible 
Navier-Stokes equations for the color-blind fluid. In addition, our new model makes 
provision for variable densities and viscosities by use of the freedom of the rest parti- 
cle equilibrium distribution, and a space dependent relaxation process. For simplicity 
we present a two-phase, two-dimensional fluid model on a hexagonal lattice in this 
paper. The extension of the current model to multi-phase fluids and other lattices, 
including a two-dimensional square lattice, a three-dimensional Face-Centered Hy- 



percubic lattice^ and a Body-Centered-Cubic JLj lattice, will be presented in detail 
in another paper |13||. 



2 Numerical Model 

Denote /,(x, t), (x, t) and as the particle distribution functions at space x and 
time t for total, red, and blue fluids respectively. Here i = 0, 1, • • •, N, where N is the 



number of moving particle directions (6 for the hexagonal lattice and 8 for the square 
lattice in two dimensions), and /j = + fj; . The lattice Boltzmann equation for 
both red and blue fluids can be written as follows: 

^(x + e il t+l) = /| fc (x,t)+Of(x,t) l (1) 

were k denotes either the red or blue fluid, and ft* = (fif) 1 + (^f) 2 is the collision 
operator. Note that in the two-phase lattice Boltzmann model by Gunstensen et al 
only a color-blind kinetic equation was given. The first term of the collision operator, 
(fi^) 1 , represents the process of relaxation to local equilibrium. For simplicity, we use 



a linearized collision operator with a single time relaxation parameter TkvA 



(n*)i = I(/*_/*W). 

Here /* is the local equilibrium state depending on the local density and velocity, 
and Tfc is the characteristic relaxation time for species k. In lattice gas automaton 
models, the form of the local equilibrium state and rate of relaxation to the local 
equilibrium distribution are completely determined by the particle scattering process. 
From a statistical point of view, the macroscopic effects of these collisions only arise 
in the transport coefficients. Therefore, lattice Boltzmann models ignore particle 
collision details such as collision cross sections and frequencies. In addition, the local 



equilibrium state can be arbitrarily chosen |TT|, with the exception that it must satisfy 
the conservation of mass and momentum: 

i i 



b(eq) 
i > 



and 



k(eq) 
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Here p r and pb are densities of the red and blue fluids respectively, p = p r + pb is the 
total density and v is the local velocity. 

Following the pressure-corrected Galilean-invariant single-phase lattice Boltzmann 



model proposed by reference |1.TJ, we assume the following equilibrium distribution 
for both red and blue fluids: 

fi (eq) = p r /(6 + m r ) + p r (e i -v/3 + 2/3(e i ) a (e i )f 3 v a Vf 3 -^v 2 ), 



fo {eq) = p r m r / (6 + m r ) - p r v 2 , 

fi (eq) = p b /(G + m b ) + pb(e i -\/3 + 2/3(e i ) a (e i ) f3 v a v f3 -^-v 2 ), 

fo (eg) = p b m b /(6 + m b )-pbV 2 . (2) 

Here we have introduced two parameters, m r and m&. They can be understood 
as the number of degenerate rest particle states for the red and blue fluids re- 
spectively To achieve a stable interface, we furthermore assume that the moving 
particle distributions for both red and blue fluids are the same when v = 0:, i.e., 
d = p r /(6 + m r ) = p b / (6 + m b ). This implies the following density ratio: ^ = fq^- 
The second part of the collision operator is similar to that given in Gunstensen et 
a/'s paper: 

(^) 2 = ^|F|(( ej -F) 2 /|F| 2 -l/2), 
where F is the local color gradient, defined as: 

F(x) = J2 e i(Pr(x + a) - p fe (x + a)). 

i 

Note that in a single-phase region of our incompressible fluid model, F = 0. Thus 
the second term of the collision operator, (Qf ) 2 , only has contribution at two-phase 
interfaces. The parameter Ak is a free parameter, controlling the surface tension 
(shown before). For simplicity, we choose A r = Ab in this paper. 

To obtain a surface tension, we follow Rothman's scheme Ol to redistribute colored 



particles at two-phase interfaces (without changing the total particle distribution, fc) 
by enforcing the red color momentum j r = J2i f[ e i, to align with the direction of the 
local color gradient. In other words, we redistribute the red density at an interface 
to maximize the following quantity: 

-CT-F). 

The blue particle distribution can then be recovered using: j\ — fa — f[. 

To derive hydrodynamics, we use a long-wavelength, low-frequency approximation 
and a multiscaling analysis as follows: 

d d 2 d 



dt dh dt 

d d 
e- 
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dx dxi 

Here t\ and t 2 represent fast and slow time scales respectively, and e is assumed to 
be a small expansion parameter. 

A Taylor series expansion of Equation (]!]) to second order in the lattice spacing 
and time step gives the following continuum kinetic equation: 

M + „ . v/ f + ie ( e< : VV/f + e, . v|/f + i||yf = Of- (3) 

Taking the zero and first order moments of over equation (^j, we readily obtain 
the continuum equations for the red and blue fluids, and the following momentum 
equation for the color-blind (total) fluid: 

^ + £e«e,(/f<-> + r>) = 0, 
where is the next order perturbation for f^. Note that since 

E(^ fc ) 2 = £(^)% = o 
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and 



F = ^V(p r -p 6 ) + e + o(e 2 ), 



the second term in the collision operator does not contribute to the continuum and 
momentum equations of the first order approximation. It, however, will contribute 
to the pressure term at an interface: P = Pq + e|VF| , where Pq comes from the 
equation of state. 

Our two-phase immiscible fluid model contains three fluid regions: a red or blue 
homogeneous region and a thin region where the two fluids mix at the interfaces. 
In a homogeneous region, the evolution of our model will recover the Navier-Stokes 
equations with viscosity of u k = (2r fe — l)/8 and sound speed of = ^3/(6 + rrik) 



according to the Chapman- Enskog expansion shown in reference ||12|| . A variation of 
viscosities for the two fluids can be obtained by choosing different r^. 

The last issue remaining to be addressed is the interfacial dynamics that take place 
in a region where red and blue fluids are adjacent. Usually the thickness of such an 
interface will depend on an averaged relaxation time, 2r r Tb/(r r + rj), and the rest 
particle distribution. There are several ways to construct a relaxation parameter so 
that the red and blue fluids have a smooth change of viscosity at their interfaces. 
To do this, we define an order parameter, if), depending on red and blue densities as 
follows: ip = p p r + p p b b • Note that in general < 1, however in a purely red fluid region, 
■0 = 1, and in a purely blue fluid region, ip = — 1. To continuously connect relaxation 
parameters r r and r& at an interface, we employ the following simple formula: 

T r ij) > 8 

g r (iff) 5>ip>0 
gbty) > > -5 
n ip < -8 

where g r (ip) and gb{ip) are second order functions of i/j: 

g r (ijj) = a + f3ijj + 7^ 2 , 
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and 

g b (i[>) = a + rjip + ^ip 2 . 

Assuming that g T {8) = T r , <5) — r&, |^ = at |^| = 5 and that g r (0) = g b (0), we 
have: a = 2r r T b / (r r + n), (3 = 2(a-r b )/8, 7 = (3/(28), r] = 2(r r -oc)/8 and £ = ^/(2<5). 
Here 5 is a free parameter controlling the interface thickness, 8 < 1. 

Using the standard definition of surface tension 0, a = J(P n — Pt)dx, and after 
some algebra jT3|, we obtain a theoretical formula for the surface tension on a two- 
dimensional hexagonal lattice: a = 9 A < r > d(12 + m r + m b ), where < r > is 
the averaged relaxation time step across the interface. Here the integration is over 
whole space along the direction perpendicular to a given interface, and P n and P t 
are respectively the normal and tangential stress tensor at the two-phase interface. 
For the interface relaxation interpolation scheme in this paper, we use < r >= 
2r r rfe/(r r + r b ). In Fig. 1, we show the theoretical prediction ( — ) and numerical 
measurements (o and +) of a as a function of particle density d with a mass ratio, 
m r /m b = 2. The o symbols represent numerical measurements of surface tension 
using the above definition, and the + signs represent the surface tension obtained by 
Laplace's formula with measurements of the pressure drop across a bubble. As seen, 
the theoretical prediction and numerical measurements agree very well in both cases. 

3 Numerical Simulations and Discussion 

To demonstrate the application of the presented lattice Boltzmann model, we show 
three numerical examples in the following sections. These simulations have been 
carried out using the CM-200 at the Advanced Computing Laboratory at Los Alamos 
National Laboratory. The first numerical test is the study of the phase segregation 
process by two-dimensional spinodal decomposition when the fluids have differing 
densities. Fig. 2 is a snapshot (t = 8300) of the two-phase area distribution when 
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the density ratio, 7 = is 10. One can see that the red fluid represents a high 
density region and the blue fluid represents a low density region. The initial condition 
for this simulation is constant density and random color distribution with periodic 
boundaries. Here the surface tension parameter is set to A = 0.01. The current model 
preserves the basic two-phase segregation processes and, in simulations not shown 
here, we have seen the preservation of stable two-phase interfaces for much higher 
density ratios (~ 200). This capability will allow us to approach the simulation of 
gas and oil mixing flows. 

The second numerical test is two-dimensional channel flow with two fluids having 
different viscosities. A lattice length of 65 in the y-direction (W = 65) and 128 in 
the x-direction (L = 128) is used, with a non-slip condition on both the upper and 
bottom boundaries. Initially, the lower half space (from lattice 1 to 32) is filled with 
red fluid, and the top half space (from lattice 34 to 65) is filled with blue fluid. The 
middle line is initialized to half red and half blue, and the relaxation parameters for 
red and blue fluids are assigned values of r r = 2 and r& = 1. Using the viscosity 
formula discussed above, this translates into a kinematic viscosity ratio of 3:1 for the 



is 



red and blue fluids respectively. A forcing technique, as described in reference |12 
used to establish a small pressure gradient across the length of the channel. Assuming 
a thin interface in the middle, an analytical solution of the velocity distribution as a 
function of y can be derived, centered on the middle of the channel fli~5"H : 



v 



(Po - P L )w 2 
2/xiL 
b _ (Pp - P L )w 2 
2/i 2 L 



2fi r u r - fi b y y 



2 fi b + /i r - Hby y 2 



(4) 

where (P ~ Pl)/L is the pressure gradient across the channel, \i T and \i b are shear 
viscosities for red (0.7875) and blue (0.2625) fluids respectively. Here also w = \^3W/4 
is half the channel width, where the factor V3/2 is a hexagonal lattice effect. In Fig. 
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3 we present the analytical prediction ( — ) and direct numerical simulation (o) for 
velocity as a function of channel width y. As seen, the agreement is excellent. 

The last numerical test is a two dimensional Hele-Shaw viscous fingering experi- 
ment. The upper and lower walls are assigned a no-slip boundary condition, and the 
fluids are given a viscosity ratio of 1:10 by assigning r r = 1 and = 5.5. To develop 
a Hele-Shaw pattern, an initial perturbation is introduced at the interface |T6|, [T7] 
which is then forced by maintaining a high pressure at the inlet and low pressure 
at the outlet. The shape of the finger can be altered by a change in the surface 
tension parameter A, given as 0.0065 in this experiment. Fig. 4 shows the evolution 
of a less viscous fluid penetrating one having a higher viscosity at times = 0, 2000, 
4000, 6000, 8000 and 10000. We see that a stable finger develops and is maintained 
throughout the simulation, although we expect a tip-splitting instability to occur for 
smaller values of A, as demonstrated by reference ||17|| . This phenomenon will be 
discussed in more detail in a subsequent paper. The fluid patterns shown here are 



similar to results obtained by other numerical methods|T6|, F7|. 

The lattice Boltzmann model in this paper is the extension of the model proposed 
by Gunstensen et al. The current model gives an exact Navier-Stokes solution, in the 
incompressible limit, for each fluid individually. The surface tension at an interface 
satisfies the Laplace formula, and the model has the capability to simulate two-phase 
fluid flows with different densities and viscosities. The simulation results for several 
typical two-phase fluid flows are shown, in good agreement with analytical solutions 
and some experimental observations. Several issues, however, remain to be addressed. 
First, to achieve a density variation, the current model uses the freedom of the rest 
particle equilibrium distribution. Although the scheme is simple, it is difficult to 
obtain high velocity flows since the assignment of a high particle mass ratio will 
accumulate a large number of particles in the rest state, reducing the fluid speed. In 
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addition, the equation of state of the current model is that of an ideal gas similar 
to the original FHP model. To simulate fluid flows with a large density variation, a 
lattice gas model with an equation of state for non-ideal gases (liquids) may provide 
some alternatives. Second, the current model uses two parabolic curves to match the 
interface relaxation characteristic time in the mixing region. For a very thin interface 
(for example, one lattice width), the parameter, 5, will be negligible. On the other 
hand, for flows with slow relaxation, the interface will be much wider. The relaxation 
form and the choice of 5 may affect interfacial dynamics. Other functional forms, 
such as a bi-normal distribution, may give better results. 

We thank G. D. Doolen and Hudong Chen for useful suggestions. Eugene Loh has 
contributed to the original CM-2 code development. This work is supported by the 
US Department of Energy at Los Alamos National Laboratory. Numerical simula- 
tion were carried out using the computational resources at the Advanced Computing 
Laboratory at the Los Alamos National Laboratory. 
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4 Figure Captions 

Fig. 1 The theoretical prediction ( — ) and numerical measurements (o and +) of 
surface tension a as a function of particle density d with a mass ratio 7 = 2. 

Fig. 2 A snapshot at time t = 8300 of the area density distribution of a spinodal 
decomposition simulation where the two fluids have a 10:1 mass ratio, 7. The 
boundaries are periodic, and the initial condition is assigned a constant density 
with a random color distribution. The surface tension parameter is assigned to 
A = 0.01. 

Fig. 3 The analytical prediction ( — ) and direct numerical simulation (o) of velocity 
as a function of channel width y. 

Fig. 4 The formation of a stable finger in a Hele-Shaw cell at times t = 0, 2000, 4000, 
6000, 8000 and 10000. A viscosity ratio of 1:10 is established by setting r r = 1 
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and r b = 5.5. The low-viscosity fluid (red) is penetrating the higher viscosity 
region (blue). 
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